#ifndef TtSemiLepSignalSel_h
#define TtSemiLepSignalSel_h


//#include "Sphericity.h"

#include "EventShapeVariables.h"


#include "TObject.h"
#include "TNamed.h"

#include "TString.h"
#include "TList.h"
#include <vector>
#include "TMath.h"
#include "TLorentzVector.h" // to perform boost

#include "MathRoot.h"

class TtSemiLepSignalSel :public TNamed 
{
  // common calculator class for likelihood
  // variables in semi leptonic ttbar decays
public:

  TtSemiLepSignalSel():TNamed("TtSemiLepSignalSel",(const char *)("")) { }
  TtSemiLepSignalSel(TString name, TString header,const std::vector<math::XYZTLorentzVector>&);

  ~TtSemiLepSignalSel();
  
  double sumEt() const { return var_sumEt; }
  double Et1() const { if (var_sumEt>0) return var_Et1/var_sumEt; return -1e10; }
  double lepeta() const { return fabs(var_lepeta); }
  double MET() const { return var_MET; }
    
  double dphiMETlepton() const { return var_dphiMETlepton; }
  
  double detajet2jet3() const { return var_detajet2jet3; }
  double detajet3jet4() const { return var_detajet3jet4; }

  double mindijetmass() const { return var_mindijetmass/massalljets; }
  double maxdijetmass() const { return var_maxdijetmass/massalljets; }

  double mindRjetlepton() const { return var_mindRjetlepton; }
 
   
  double DeltaPhi(math::XYZTLorentzVector, math::XYZTLorentzVector);
  double DeltaR(math::XYZTLorentzVector, math::XYZTLorentzVector);

double sphericity() const { return var_sphericity;}
double aplanarity() const { return var_aplanarity;}

double Et2 () const { if (var_sumEt>0) return var_Et2/var_sumEt; return -1e10;}
double Et3 () const { if (var_sumEt>0) return var_Et3/var_sumEt; return -1e10;}
double Eta1() const {return var_Eta1;}
double Eta2() const {return var_Eta2;}
double Eta3() const {return var_Eta3;}
double dptjet1jet2 () const  {if (var_sumEt>0) return var_dptjet1jet2/var_sumEt; return -1e10;} 
double dptjet1jet3 () const  {if (var_sumEt>0) return var_dptjet1jet3/var_sumEt; return -1e10;} 
double djet1jet2pt () const  {if (var_sumEt>0) return var_djet1jet2pt/var_sumEt; return -1e10;}


double detajet1jet2() const { return var_detajet1jet2; }
double dthetajet1jet2_boost() const { return var_dthetajet1jet2_boost;}
double dphijet1jet2_boost() const { return var_dphijet1jet2_boost;}

double Et2byEt1 () const {if (var_Et1>0 && var_Et2>0 ) return var_Et2/var_Et1; return -1e10;}
double Et3byEt1 () const {if (var_Et1>0 && var_Et3>0) return var_Et3/var_Et1; return -1e10;}
double Et3byEt2 () const {if (var_Et2>0 && var_Et3>0)  return var_Et3/var_Et2; return -1e10;}

double dphijet1jet2 () const { return var_dphijet1jet2;}
double dphijet1jet3 () const { return var_dphijet1jet3;}
double dphijet2jet3 () const { return var_dphijet2jet3;}

double sphericity_boost() const {return var_sphericity_boost;}

private:

double var_dphijet1jet2;
double var_dphijet2jet3;
double var_dphijet1jet3;

double var_sphericity_boost;

double var_Et2;
double var_Et3;
double var_Eta1;
double var_Eta2;
double var_Eta3;
double var_dptjet1jet2;
double var_dptjet1jet3;
double var_djet1jet2pt;
double var_detajet1jet2;
double var_dthetajet1jet2_boost;
double var_dphijet1jet2_boost;

  double var_sphericity;
  double var_aplanarity;

  
  double var_sumEt;
  double var_Et1;
  double var_lepeta;
  double var_MET;
  
  double var_dphiMETlepton; 
  
  double var_detajet2jet3;
  double var_detajet3jet4;

  double var_mindijetmass;
  double var_maxdijetmass;

  double var_mindRjetlepton;
  

  double massalljets;

public:
 	ClassDef(TtSemiLepSignalSel,2);
};


 #ifdef __MAKECINT__
ClassImp(TtSemiLepSignalSel)
#endif


#if defined(TtSemiLepSignalSel_CC)


TtSemiLepSignalSel::TtSemiLepSignalSel(TString name, TString header,const std::vector<math::XYZTLorentzVector>& topJets):TNamed(name,header)
{ 


  unsigned int nJetsMax = topJets.size();
  
  var_MET = -1e+10;
  var_sumEt = 0.;

//std::cout<<"TtSemiLepSignalSel!!! "<< var_MET<<std::endl;  

//std::cout<<nJetsMax<<std::endl;

  math::XYZTLorentzVector Jetsum(0.,0.,0.,0.);
  
  for(unsigned int i=0; i<nJetsMax; i++) {
    math::XYZTLorentzVector aJet = topJets[i];
    Jetsum += aJet;
    var_sumEt += topJets[i].Et();
  }
  massalljets = Jetsum.M();

//std::cout<<"TtSemiLepSignalSel!!! "<<massalljets<<std::endl;  

 
  var_lepeta = -1e+10;

//std::cout<<"TtSemiLepSignalSel!!! "<<var_lepeta<<std::endl;  


//  math::XYZTLorentzVector Met = MET->p4();
  std::vector<double> Etjet;
  std::vector<double> Jetjet;
  double dijetmass;
  var_mindijetmass = 99999.;
  var_maxdijetmass = -1.;

  for(unsigned int i=0; i<nJetsMax; i++) {
    math::XYZTLorentzVector aJet = topJets[i];
    Etjet.push_back( aJet.Et());
    for(unsigned int j=i+1; j<nJetsMax; j++) {
      math::XYZTLorentzVector asecJet = topJets[j];
      dijetmass = (aJet+asecJet).M();
      if(dijetmass<var_mindijetmass) var_mindijetmass = dijetmass;
      if(dijetmass>var_maxdijetmass) var_maxdijetmass = dijetmass;
    }
  }

//std::cout<<"var_mindijetmass = "<<var_mindijetmass<<std::endl;


if (Etjet.size()>0)  var_Et1 = Etjet[0];
   else var_Et1=0;
  var_dphiMETlepton =-1e+10;


std::vector<math::XYZTLorentzVector> _vec1;
std::vector<math::XYZTLorentzVector> _vec2;

for(unsigned int i=0; i<nJetsMax; i++ ) _vec1.push_back(topJets[i]);
for(unsigned int i=0; i<nJetsMax; i++ ) _vec2.push_back(topJets[i]);

if (nJetsMax>2)  var_detajet2jet3 = fabs(_vec1[1].Eta()-_vec2[2].Eta());
if (nJetsMax>3)  var_detajet3jet4 = fabs(_vec1[2].Eta()-_vec2[3].Eta());


//std::cout<<"var_detajet2jet3 = "<<var_detajet2jet3<<std::endl;

 
  var_mindRjetlepton = -1e+10;


///Event shape variables

std::vector<math::XYZVector> _vecXYZ;
for ( unsigned int i=0;i<nJetsMax;i++) _vecXYZ.push_back(topJets[i].Vect());
EventShapeVariables _etsh(_vecXYZ);

///New code! to test
///sphericity, sphericity_boost
/**
	TObjArray *  _jets = new TObjArray();
	_jets->SetOwner(kTRUE); // To 'autodelete' all daughters

	for ( unsigned int i=0;i<nJetsMax;i++) _jets->Add(new TLorentzVector(topJets[i].X(),topJets[i].Y(), topJets[i].Z(), topJets[i].T() ));
        Sphericity * spher = new Sphericity();
        spher->analyze(_jets,2);
	delete _jets;
**/		

if (topJets.size()>0) {
var_sphericity = _etsh.sphericity();
///var_sphericity = spher->sphericity();
var_aplanarity = _etsh.aplanarity();
} 
else {

var_sphericity = -1e10;
var_aplanarity = -1e10;

}
if (Etjet.size()>1)  var_Et2 = Etjet[1];
   else var_Et2=-1e10;
if (Etjet.size()>2)  var_Et3 = Etjet[2];
   else var_Et3=-1e10;

if (topJets.size()>0)  var_Eta1 = topJets[0].Eta();
   else var_Eta1=-1e10;
if (topJets.size()>1)  var_Eta2 = topJets[1].Eta();
   else var_Eta2=-1e10;
if (topJets.size()>2)  var_Eta3 = topJets[2].Eta();
   else var_Eta3=-1e10;

//std::cout<<" 1111!!!"<<std::endl;

if (topJets.size()>1)  var_dptjet1jet2 = topJets[0].Et() - topJets[1].Et();
else var_dptjet1jet2=-1e10;
if (topJets.size()>2)  var_dptjet1jet3 = topJets[0].Et() - topJets[2].Et();
else var_dptjet1jet3=-1e10;
if (topJets.size()>1)  var_djet1jet2pt = (topJets[0]- topJets[1]).Et();
else var_djet1jet2pt=-1e10;
if (nJetsMax>1)  var_detajet1jet2 = fabs(_vec1[0].Eta()-_vec2[1].Eta());
else var_detajet1jet2=-1e10;



//std::cout<<"var_sphericity = "<<var_sphericity<<std::endl;



std::vector<TLorentzVector>_vecs12;
if (topJets.size()>1) 
{
for(unsigned int i=0; i<2; i++ ) _vecs12.push_back(TLorentzVector(topJets[i].X(),topJets[i].Y(),topJets[i].Z(),topJets[i].T()));
TLorentzVector _Higgs = _vecs12[0]+_vecs12[1];
//TVector3 _CM = -_Higgs.BoostVector();
_vecs12[0].Boost( -_Higgs.BoostVector());
_vecs12[1].Boost( -_Higgs.BoostVector());

var_dthetajet1jet2_boost = ROOT::Math::VectorUtil::Angle<TLorentzVector,TLorentzVector>(_vecs12[0],_vecs12[1]); 
var_dphijet1jet2_boost = ROOT::Math::VectorUtil::DeltaPhi<TLorentzVector,TLorentzVector>(_vecs12[0],_vecs12[1]); 

//std::cout<<" _vecs12[0].Boost.X = "<< _vecs12[0].X()<<std::endl; 

} else {
var_dthetajet1jet2_boost=-1e10;
var_dphijet1jet2_boost = -1e10;


}

if (topJets.size()>1)
var_dphijet1jet2= ROOT::Math::VectorUtil::DeltaPhi<math::XYZTLorentzVector,math::XYZTLorentzVector>(topJets[0],topJets[1]);
else var_dphijet1jet2=-1e10;

if (topJets.size()>2){
var_dphijet1jet3= ROOT::Math::VectorUtil::DeltaPhi<math::XYZTLorentzVector,math::XYZTLorentzVector>(topJets[0],topJets[2]);
var_dphijet2jet3= ROOT::Math::VectorUtil::DeltaPhi<math::XYZTLorentzVector,math::XYZTLorentzVector>(topJets[1],topJets[2]);

}else {
var_dphijet1jet3=-1e10;
var_dphijet2jet3=-1e10;
}

//std::cout<<"22323 !!!"<<std::endl;

_vecs12.clear();
_vecXYZ.clear();
for(unsigned int i=0; i<nJetsMax;i++) _vecs12.push_back(TLorentzVector(topJets[i].X(),topJets[i].Y(),topJets[i].Z(),topJets[i].T()));
TLorentzVector _CM (0.,0.,0.,0.); 
for(unsigned int i=0; i<nJetsMax;i++) _CM+=_vecs12[i];

//cout<<"Before boost"<<endl;

for(unsigned int i=0; i<nJetsMax;i++) _vecs12[i].Boost(-_CM.BoostVector());


for(unsigned int i=0; i<nJetsMax;i++) _vecXYZ.push_back(math::XYZVector(_vecs12[i].X(),_vecs12[i].Y(),_vecs12[i].Z()));
EventShapeVariables _etsh_boost(_vecXYZ);

//cout<<"After boost"<<endl;

///New code! to test
///sphericity, sphericity_boost
/**
        TObjArray *  _jets2 = new TObjArray();
        _jets2->SetOwner(kTRUE); // To 'autodelete' all daughters
        for ( unsigned int i=0;i<nJetsMax;i++) _jets2->Add(new TLorentzVector(_vecs12[i].X(),_vecs12[i].Y(), _vecs12[i].Z(), _vecs12[i].T() ));
  
        spher->analyze(_jets2,2);
	delete _jets2;
**/


if (topJets.size()>0)
var_sphericity_boost = _etsh_boost.sphericity();
///var_sphericity_boost = spher->sphericity();
else
var_sphericity_boost = -1e10;

//std::cout<<"END "<<std::endl;

///if (spher) delete spher;


  
}




double TtSemiLepSignalSel::DeltaPhi(math::XYZTLorentzVector v1, math::XYZTLorentzVector v2)
{
  double dPhi = fabs(v1.Phi() - v2.Phi());
  if (dPhi > TMath::Pi()) dPhi =  2*TMath::Pi() - dPhi;
  return dPhi;
}

double TtSemiLepSignalSel::DeltaR(math::XYZTLorentzVector v1, math::XYZTLorentzVector v2)
{
  double dPhi = DeltaPhi(v1,v2);
  double dR = TMath::Sqrt((v1.Eta()-v2.Eta())*(v1.Eta()-v2.Eta())+dPhi*dPhi);
  return dR;
}

TtSemiLepSignalSel::~TtSemiLepSignalSel() 
{
}


#endif

#endif

